#Alexander F. Gazmararian
#afg2@princeton.edu
#January 9, 2024

#Code to create Figure 1:
#\textbf{Shale Shock Began in 2008.}
#The left panel shows normalized coal and gas production in the US.
#The right panel shows net electricity generation.
#Source: EIA

#Load packages
library(tidyverse)
library(tidylog)
library(gridExtra)
library(here)

#Load data
##Run: "01_clean_naturalgasannual.R" to create this dataset
ng <- readRDS(here("data", "input", "eia_naturalgasannual/eia_nga.rds"))
##Run: "01_clean_eia_electricityfuel.R" to create this dataset
gen <- readRDS(here("data", "input", "eia_monthlyenergyreview", "eia_mer.rds"))
##Run: "01_clean_eia_annualcoalreport.R" to create this dataset
coal <- readRDS(here("data", "input", "eia_annualcoalreport", "eia_annualcoalreport.rds"))

#Merge natural gas and coal production data into same frame
prod <- full_join(ng, coal, by = "year")
#create total coal variable
prod$coal <- with(prod, underground + surface)
#Create function for min-max normalization
min_max_norm <- function(x) {(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))}
#Normalize the variables
prod$coal_n <- min_max_norm(prod$coal)
summary(prod$coal_n)
prod$gas_n <- min_max_norm(prod$gasprod)
#Create plot of normalized production
plot.prod <- prod %>%
  ggplot() +
  geom_vline(xintercept = 2008, color = "red", size = 0.5) +
  geom_line(aes(x = year, y = gas_n), color = "gray", size = 1.25, lty = "dashed") +
  geom_line(aes(x = year, y = coal_n), color = "black", size = 1.25) +
  scale_x_continuous(breaks=seq(1940,2020,10)) +
  scale_y_continuous(limits = c(0, 1.04)) +
  labs(x = "Year", y = "Production (normalized)", title = "Production") +
  annotate("text", x = 2020, y = .1, label = "Coal", color = "black", size = 6) +
  annotate("text", x = 2020, y = 1.04, label = "Gas", color = "gray", size = 6) +
  theme_bw(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust=.5)
  )
#Create plot of net-generation
plot.gen <- gen %>%
  ggplot() +
  geom_vline(xintercept = 2008, color = "red", size = 0.5) +
  geom_line(aes(x = year, y = netgen_NaturalGas / netgen_Total), color = "gray", size = 1.25, lty = "dashed") +
  geom_line(aes(x = year, y = netgen_Coal / netgen_Total), color = "black", size = 1.25) +
  scale_x_continuous(breaks = seq(1950, 2020, 10)) +
  scale_y_continuous(limits = c(0, 1.04)) +
  labs(x = "Year", y = "Net-Generation Share", title = "Net-Generation by Fuel Source") +
  annotate("text", x = 2020, y = .15, label = "Coal", color = "black", size = 6) +
  annotate("text", x = 2020, y = .46, label = "Gas", color = "gray", size = 6) +
  theme_bw(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust=.5)
  )
#Create combined figure
plot.out <- grid.arrange(plot.prod, plot.gen, ncol = 2)
plot.out
#Save plot output
ggsave(
  device = "tiff",
  filename = here("output", "figures", "fig1_shocktiming.tiff"),
  plot.out,
  scale = 1.5,
  width = 6.5,
  height = 2.5,
  dpi = 300
)
ggsave(
  filename = here("output", "figures", "fig1_shocktiming.png"),
  plot.out,
  scale = 1.5,
  width = 6.5,
  height = 2.5,
  dpi = 300
)
